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We explore to what extent path-integral quantum Monte Carlo methods can efficiently simulate 
the tunneling behavior of quantum adiabatic optimization algorithms. Specifically we look at sym¬ 
metric cost functions defined over n bits with a single potential barrier that a successful optimization 
algorithm will have to tunnel through. The height and width of this barrier depend on n, and by 
tuning these dependencies, we can make the optimization algorithm succeed or fail in polynomial 
time. In this article we compare the strength of quantum adiabatic tunneling with that of path- 
integral quantum Monte Carlo methods. We find numerical evidence that quantum Monte Carlo 
algorithms will succeed in the same regimes where quantum adiabatic optimization succeeds. 


I. INTRODUCTION 
A. Background 

Quantum Adiabatic Optimization (qao), first pro¬ 
posed by Farhi et al. [lj, is a quantum algorithm for deter¬ 
mining the minimum of a cost function by slowly evolv¬ 
ing a Hamiltonian from one with known ground state to 
one that has as its ground state the solution to an opti¬ 
mization problem. QAO relies on the quantum adiabatic 
theorem (see Jansen et al. [z] for a proof), which roughly 
says that a system is guaranteed to stay in its ground 
state if the Hamiltonian evolution time-scales are much 
larger than the square of the inverse spectral gap. It was 
shown by, for example, Reichardt Q and Farhi et al. Q 
that this algorithm might provide an exponential speed¬ 
up over naive local search algorithms. On the other hand, 
Farhi et al. 0 also described cases where it is no better 
than the quadratic speed up of Grover’s quantum search. 

QAO is sometimes referred to as quantum annealing 
(qa) in relation to classical simulated annealing (sa) 
where a simulated system’s temperature is slowly lowered 
reducing the probability of energetically less favorable 
states until the ground state is reached at zero tempera¬ 
ture. In recent years, there has been a push to compare 
these two methods, often by using simulated quantum 
annealing (sqa) [§]. SQA uses a path integral expan¬ 
sion of the partition function for the evolving system to 
create a (d + 1) dimensional classical system on which 
Monte Carlo techniques can be used. Instead of varying 
the temperature as in SA, SQA varies the Hamiltonian in 
the same way as QAO. 

This path-integral Quantum Monte Carlo (qmc) al¬ 
gorithm has been used to compare classical SA and QA. 
Heim et al. [3] among others have shown that QMC meth¬ 
ods outperform classical SA in several cases. In other sit¬ 
uations Battaglia et al. 0 showed that SA can perform 
better than QMC. In addition, new techniques in SQA 
through QMC continue to be developed and improved, 


such as by Farhi et al. 

SQA through QMC captures much of the power of QAO, 
and for some problems these two methods show corre¬ 
lation in their success rates while at the same time re¬ 
maining uncorrelated from classical SA 0 - However, 
Hastings has recently 0 constructed several examples 
where QAO will find the ground state in polynomial time 
whereas QMC methods will take exponential time. 

B. Central Problem 

To directly compare the strengths and weaknesses QAO 
and QMC methods, this article will look at their respective 
efficiencies as the two methods tunnel through a potential 
barrier. The specific problem consists of a symmetric cost 
function on n bits where each basis state x £ {0, l} n is 
weighted by its Hamming weight |x| in combination with 
a potential barrier centered at \x\ = n/ 4. Barriers of 
this form have been partially considered in the context of 
QAO by Reichardt Q who found that QAO would succeed 
in time polynomial in n if the height and width of the 
barrier are both ^(n 1 / 4 ). 

A simplified problem with a barrier of width 1 was 
analyzed by Farhi et al. 3, comparing QAO and classi¬ 
cal SA. There QAO was found to succeed in polynomial 
time while classical SA could not. Crosson and Deng jl2| 
showed that the same thin barrier limit is a case where 
QMC methods and QAO both succeed together. 

Muthukrishnan et al. [l3| analyzed a similar problem 
where instead of a barrier in the Hamming Weight, they 
have a plateau. This problem has a constant gap, but 
they showed that QAO still outperforms SA, though both 
run in polynomial time. They also showed that a non- 
adiabatic approach to QA could outperform QAO in this 
case with a constant gap. 

Our current goal is to extend the comparison of QAO 
and QMC methods to the case of a varying barrier size. 
Therefore, we seek to determine if the correlation be¬ 
tween the two continues for the full case where the width 


2 


and height of the barrier are both powers of the number 
of qubits n. 

C. Organization 


where i sums over all n qubits. The ground state of this 
Hamiltonian is a uniform superposition over all |x) states. 
Therefore, the ground state is initially a binomial proba¬ 
bility distribution over \x\ with width ~ yfn centered at 
|x| = n/2. In QAO, we create the Hamiltonian 


In Section [TTJ we will setup the particular problem we 
are working with, defining our symmetric Hamiltonian 
and its tunable parameters. In Section 1IIII we will ex¬ 
amine the energy eigenvalues of this Hamiltonian. We 
will focus on the spectral gap between the ground state 
and first excited state and will primarily use numerical 
diagonalization. The size of this spectral gap determines 
how slowly adiabatic evolution must go in order to stay 
in the ground state. 

Section llVl will outline and develop on our Monte Carlo 
method. We will go through the approximations and how 
those approximations effect our simulations; additionally, 
we will discuss our choice of update rules. In Section lYl 
we present the results of our Monte Carlo simulations 
and compare the scaling behavior of these simulations to 
the scaling behavior of the spectral gap from Section Hill 
Finally in Section IVII we discuss the limitations of our 
Quantum Monte Carlo algorithm and present several av¬ 
enues for extension and generalization of our work. 


II. HAMMING WEIGHT WITH A BARRIER 


Our problem is one discussed by Reichardt [ 3 }, and a 
simplified version of it was analyzed by Crosson and Deng 
[T3 . We consider a symmetric cost function /(|x|) = 
|x|+6(|x|), where |x| is the Hamming Weight of the length 
n bit string x, and b(z) is some perturbing function. We 
will take b(z ) to be some barrier, centered around z = 
n/4, that has width and height proportional to n a . For 
ease of computation, we will use 


b(z) 


n a when (f - ±cn Q ) < z < (f + \cn a ) 
0 otherwise 


(1) 

where c is an n independent constant. From now on we 
will say that this barrier has size cn a . The full cost 
function will have a global minimum at |x| = 0 and a 
local minimum at \x\ = \_j + \cn a \ + 1. 

We will encode this problem into a Hamiltonian on a 
Hilbert space of n qubits: 


H i= E /(\ x \)\ x )( x \- ( 2 ) 

xG{0,l} n 


In QAO, we slowly transition from a Hamiltonian with 
a known ground state into one with a desired ground 
state such as Hi (e.g. in this problem, we want to find 
the \x\ = 0 state). The standard initial Hamiltonian is 


H 0 = £(Ho)i 
2—1 


with 



( 3 ) 


H(s) = (l-s)H 0 + sH 1 , (4) 

where s goes from 0 to 1. If we vary s slowly enough, 
the adiabatic theorem says that the system will remain 
in the ground state. Therefore, the system will be forced 
to tunnel through the potential barrier in order to reach 
the true final ground state with |x| = 0. As s changes, 
the first two energy eigenlevels remain distinct and have 
some spectral gap g(s). If the minimum gap over s is 
min sg r 0 .i] g(s) = g m in, then adiabatic evolution is guar¬ 
anteed to keep the system in the ground state if it takes 
time 


III. EXACT SPECTRAL GAP 

To determine the minimum spectral gap g m i n we will 
numerically diagonalize the Hamiltonians H(s). Using 
the symmetry of the Hamiltonians we are able to do this 
accurately in the same range of finite n that our Quantum 
Monte Carlo simulations access. As a result, we will be 
able to compare the QMC run-times directly to the 1 / g? nm 
quantity, rather than having to rely on extrapolations to 
large n behavior. 


A. Symmetrized Hamiltonian 


In order to diagonalize the H(s ) of Eq. [I] for sizable n, 
we rely on the symmetric subspace of our system. For 
each Hamming weight 0 < h < n we have that H(s) is 
degenerate in the (?) dimensional subspace spanned by 
the vectors {|x) : |x| = h}. Hence we see that the spec¬ 
trum of H(s) has at most n+1 distinct eigenvalues, which 
will simplify our numerical calculations significantly. We 
rewrite the Hamiltonian as follows. 


Hsym(s) — y] 
h —0 


+E 


h =0 


+E 


h—0 




s/O) 


0)01 


- ( \^V(h + l)(n-h) 

2 ^V o+oo-fe) 


( 5 ) 

0)0 + 1 | 

0 + 1)01 


The spectral gap is then found by diagonalizing the re¬ 
sulting (n+ 1) x (n+ 1) tridiagonal matrix. Incidentally, 
the Hilbert space of this symmetrized system is identi¬ 
cal to that of a single spin-n/2 particle. In that spin 
context, the Hamiltonian describes applying a magnetic 
field to the particle where the field starts as a uniform 
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FIG. 1: (/min vs. n for barrier size n 0 ’ 5 : We show a best 
fit linear regression through the log-log data and plot the 
residuals of that linear fit versus the log-log data. The fact 
that the residuals curve down means that </ m i n is decreasing 
faster than a power law with n, indicating superpolynomial 
growth in the QAO run-time 
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Number of Qubits 

FIG. 2: (/ m in vs. n for barrier size n 0 ' 4 : The best fit linear 
regression to the log-log data has residuals that curve down¬ 
wards indicating superpolynomial growth in the QAO run¬ 
time. Also, notice that y -axis scale on the residual plot is 
much smaller than in Fig. [l] indicating that this scaling is 
not as strong as n the higher a case. 


field in the —x direction and then rotates to one in the 
z direction with certain momentum modes picked out as 
more energetic. 

Furthermore, the adiabatic theorem states that QAO 
run-time depends on the minimum spectral gap as s 
evolves, g m i n , so we minimize the gap as a function of 
s. We restrict ourselves to n divisible by 4 so that the 
barrier is centered on an integer Hamming weight. Since 
the barrier width increases in integer steps, we only con¬ 
sider ns such that the width has just increased (i.e. n 
such that [1 + cn Q J > [1 + c(n — 4) Q J). Finally we re¬ 
quire the barrier to have width less than n/2 , preferably 
much less, so that the s = 0 ground state does not have 
a significant overlap with the region of the barrier. 


B. Numerical Results 

In Fig. [T]we show the minimum gap for a barrier of size 
n 0 5 as a function of n. The line drawn through the points 
is a linear best fit to the log-log data, and the plot below 
shows the residuals for this fit. Since the residuals curve 
downward, the gap is decreasing faster than a power law 
can account for; therefore, the running time for QAO, 
which depends on the gap <7~f n , is superpolynomial in n 
for a = 0.5. In Figs. [2] and [3J we show similar plots for 
a = 0.4 and 0.3 respectively. 

Varying a, we do the same procedure, sweeping 
through a range of n from 100 to 5000 with c = 1. The 
second derivative of these log-log residuals can be used 
to estimate the curvature of those residuals (i.e. whether 
they are concave up or down), and these second deriva¬ 
tives are plotted in Fig. [I] The residual plots in Figs. 03 [2j 
and[3]are all used in the construction of Fig. [31 Since the 
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FIG. 3: (/min vs. 71 for barrier size n 0 ' 3 : The best fit linear 
regression to the log-log data has residuals that are essentially 
zero, indicating polynomial scaling with n. We have used the 
same residual scale as in Fig. [2] to indicate just how small 
these residuals are. 
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second derivative varies over the range of n, the second 
derivative is averaged and the standard error is used as 
the error bars. A negative number indicates superpoly¬ 
nomial running time, whereas zero represents polynomial 
scaling. 

The curvature in Fig.[4]becomes negative by more than 
one error bar starting at a = 0.34, which indicates that 
the quantum adiabatic algorithm will undergo a transi¬ 
tion from polynomial to exponential scaling somewhere 
between a = 0.33 and 0.34. 

It is a folklore result 0 that as n grows for a barrier 
with height and width proportional to n a , the spectral 
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FIG. 4: Deviation of g m i n from Power Law in n: For each 
barrier scaling power a at c = 1 we found the spectral gap 
for n between 100 and 5000 and tried to fit a linear curve to 
the log-log plot of spectral gap versus n. What is displayed 
here is the curvature of the residuals from those fits. If the 
residuals are concave down (meaning negative curvature on 
this figure), the spectral gap is decreasing faster than a power 
law in n. Therefore, QAO will become superpolynomial in n 
somewhere between a = 0.33 and 0.34. 
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FIG. 5: (; m in vs. n for barrier size 2n 0 ' 4 : The best fit linear 
regression to the log-log data has residuals that curve upwards 
indicating polynomial or subpolynomial decrease with n. At 
the end of the n range, the residuals begin to curve down 
again, indicating the beginning of the superpolynomial region 
indicated by Figs. [|] and 0 


IV. PATH-INTEGRAL QUANTUM MONTE 
CARLO 


gap decreases asymptotically as: 


j 

f constant 

if a < \ 


Qmin — > 

. l/polynomial(n) 

1 1 /exponential(n) 

if \ < a < | 
if i < a 

(6) 

Hence we expect a transition from polynomial to 
nentially small gaps to occur when a = 1/3. This 

expo- 

result 


meshes exactly with the results of our simulations. Our 
numerical results are still useful in their own rights since 
our QMC calculations will be accessing finite n values and 
it is important to compare the QMC results with equiv¬ 
alent gap results, and we need these results to be aware 
of any possible small n phenomena. 

Additional numerical results indicate that the large n 
scaling behavior in Fig. 3] does not hold for smaller n 
when c is large. For instance Fig. [2] does display the 
large n superpolynomial behavior, but Fig. [5] does not. 
In Fig. 0 if we consider just the largest n, there are 
indications that the residuals are becoming concave down 
at the end, indicating that the superpolynomial scaling 
is starting at the end of the n range we are looking at. 

The computational limits of our QMC algorithm and 
computing facility mean that some of the QMC simula¬ 
tions in this article will be at lower n where the large n 
scaling behavior is not yet dominant. In cases where we 
can access the large n scaling behavior, such as a = 0.4 
and c = 1 in Fig. [2] we will mention so in subsequent 
analysis. Largely, we will be comparing QMC running 
times with g“ in directly so that we can see if QMC run¬ 
ning time scales polynomially with QAO running time. 


The path-integral QMC algorithm [l5| is a method of 
simulating a quantum mechanical system at finite inverse 
temperature /3. The procedure uses Trotter expansion to 
take an n qubit quantum system to a classical system of 
n bits evolving in a discretized “imaginary time” dimen¬ 
sion. These time evolving states can then be treated as 
states in a Monte Carlo simulation that samples possible 
paths of the system. 

The Monte Carlo algorithm then picks paths with 
probability proportional to their Boltzmann weights, so 
from these states, an expectation value for the ground 
state energy can be obtained. We run the Monte Carlo 
algorithm for fixed s until we reach the ground state at 
that s value and then transition to a new s. This so 
called annealing schedule captures the same adiabaticity 
that makes QAO so powerful. 


A. Trotter Expansion 


To start, we take the partition function at finite inverse 
temperature /3 and Trotter expand it into T “time”-slices 


Z = Tr 

E 


(7) 


= lim 

T —too 


c<°). 




Y[{x^\e~^ H \x {T+1) 

,T = 0 


where the sums go over each x^ £ {0,1}". In order 
to be in the ground state, the temperature needs to be 
low, which means high /?, but T also needs to be much 
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greater than /3 in order for the Trotter approximation to 
work well. In practice, we will take /3 = 32 and T oc n 
for reasons that will be discussed in subsection IIV Bl We 
also have periodic boundary conditions x^ = x^. The 
goal is to have the operators act on these \x ) basis states 
so that we can get a partition function in terms of c- 
numbers. Each of the T bases corresponds to a different 
imaginary “time” slice of the system, so we are transform¬ 
ing our n qubit system into an n x T lattice of classical 
bits with interactions between adjacent time slices. 


B. Exponential Approximation 

The Hamiltonian includes terms diagonal in the com¬ 
putational basis, which we will call Hd , and off-diagonal 
terms, which we will call H 0 . The goal is to separate 
out these terms so that each operator can act on its own 
eigenbasis. There are two approximations that can be 
used here: either a linear approximation or an exponen¬ 
tial approximation for fi/T —> 0: 

e -^H d+Ho ) = p _ l( Hd + Hq) + 0 (((3/T) 2 ) (8a) 

e -^(H d +H a ) _ e -±H d( ,—§,H 0 + o((/3/T) 2 ). (8b) 

To first order these are both the same, but the additional 
terms in the exponential change the algorithm signifi¬ 
cantly. The linear approximation only includes one copy 
of the off-diagonal Hamiltonian, so adjacent Trotter time 
slices would differ by at most a single bit. Single bit flips 
between adjacent sites lend a nice sense of continuity to 
the time dimension, but they necessitate larger T. The 
off-diagonal part of the Hamiltonian manifests itself in 
the simulation as bit flips between adjacent time-slices, 
so in order to get enough bit flips in the linear approxima¬ 
tion, T must be larger, whereas the exponential approx¬ 
imation, with multiple adjacent bit flips, can be more 
compact. 

In Fye [3, there is a discussion of these two approxi¬ 
mation methods where they find that for local Hamilto¬ 
nians the exponential approximation is more robust and 
can be used with an n-independent T. The linear approx¬ 
imation requires T to increase with increasing n, making 
it less desirable. Our Hamiltonian relies on the Hamming 
Weight, which is a non-local quantity, so these results do 
not hold perfectly. We found that the exponential ap¬ 
proximation did require T to have some dependence on 
n; however, numerically, we found that dependence to be 
much smaller than the dependence of the linear approx¬ 
imation. Therefore, we use the exponential approxima¬ 
tion in this article. 

Eventually, we will want to interpret the product of 
these exponentials as a Boltzmann factor or probability 
for the given n x T configuration of the system. In order 
to do this, the Boltzmann factors must be positive. In 
order to ensure that our approximated exponentials re¬ 
main positive, the Hamiltonian must be one with “no sign 


problem.” This means that all the off-diagonal terms in 
the Hamiltonian must be non-positive. To see why, con¬ 
sider Eq. [Sal if the off-diagonal Hamiltonian contained 
negative terms, then this operator would lead to nega¬ 
tive terms if it were between non-identical states. This 
same logic is true in Eq. I8bl Our Hamiltonian has no 
sign problem, so we are free to use these methods. 

C. Final Partition Function 

For an in depth derivation of the partition function see 
Appendix[Aj Here, we will just cite the resulting partition 
function 

T—1 

e -£((i-«)9+»/(l* <T) D) ( 9 ) 

r —0 
n 

/gT 2 -f (—1 ) x d X d e T 2 J . 

d—1 

The first summation can be thought of as a sum over 
possible states, where a state is a full configuration of 
the nxT bit lattice. The expression in the square brack¬ 
ets is the Boltzmann factor for that configuration. The 
Boltzmann factors are the unnormalized probabilities for 
the states, so they can be used in a Metropolis algo¬ 
rithm to create a Monte Carlo simulation. The Quan¬ 
tum Monte Carlo method consists of performing standard 
Monte Carlo methods on this classical partition function 
which can then be used to gain information about the 
original quantum system (e.g. see Appendix [B] for how 
to extract the energy from this Monte Carlo simulation). 

D. Update Rules 

We follow the same update rule as Crosson and Deng 
fl2| , where we sweep through these nxT bits. One sweep 
consists ofnxT updates, where we go through each bit 
in the lattice separately. For that bit we try flipping its 
value, and then compare the Boltzmann weight of the 
lattice before and after the bit-flip. The acceptance rate 
of this bit flip is then equal to the ratios of the Boltzmann 
factors before and after the flip. Once the sweep has gone 
through every bit in the lattice, the sweep ends, and the 
algorithm calculates the current ground state energy of 
the entire lattice based on the results of Appendix [B] 

For the annealing schedule, we have a fixed As = ^ 
and change how much time we spend on each s value. 
The algorithm calculates the quantum mechanical en¬ 
ergy (see Appendix [B]) of the system after each sweep 
and moves onto the next s value when the energy gets 
close enough to the true ground state energy. This an¬ 
nealing schedule does use information that the QMC algo¬ 
rithm would not have in a normal simulation (namely the 
ground state energy and spectral gap), but since our goal 


Z = lim 

T—>• OO 


E 


itCO). 
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FIG. 6: QMC Sweeps vs. s for barrier size 3 n 0 5 at n = 116: 
This is averaged over 30 simulations. The spike corresponds 
to tunneling through the potential barrier and is roughly at 
the location of the minimum spectral gap. 


is to judge how long it takes to reach the ground state 
rather than how long it takes the algorithm to realize it 
has reached the ground state, this is appropriate. 

The algorithm judges it is close enough to the true 
ground state when the average energy over the last 100 
sweeps, (E(s)) 100 , is within 0.4 spectral gaps, g(s), of the 
true ground state energy, Egs(s): 


\(E( s )) W0 -E GS (s) | 

9(s) 


< 0.4 


( 10 ) 


In subsequent graphs, we will report the number of 
sweeps for each s value. If the algorithm has already 
satisfied this condition after the first 100 steps, it ex¬ 
trapolates back to when it first met the update condition 
and report that as the number of sweeps. 

In Fig. [Gj the results are shown for simulations using 
a barrier of size 3n 0 ' 5 , and n = 116. Notice the spike in 
run-time corresponding to tunneling through the poten¬ 
tial barrier. The location of this spike in s corresponds 
to the location of the minimum spectral gap in the exact 
problem, and this s location always occurs in roughly the 
same location across multiple values of n and a. In the 
next section when we report the run-time of the QMC sim¬ 
ulations, we will report the total number of sweeps taken 
between s = 0.3 and s = 0.5. For all of our simulations, 
this s range captures the run-time spike and some of the 
surrounding area while ignoring any low s initialization 
artifacts or high s tailing-off. 


V. NUMERICAL MONTE CARLO RESULTS 

In this section, we will explore a few different values 
of the barrier scaling power, a, and the width scaling 
coefficient, c, using the QMC methods developed in the 
previous section. For most of the simulations considered 


FIG. 7: QMC Sweeps vs. g m in' This lists all our data together; 
a further breakdown of this data is available in Figs. [8] [9] and 
m Notice that there is an obvious strong correlation between 
required and sufficient QMC sweeps and the gap. More analy¬ 
sis, specific to the different a values can be found in Figs. [8] 

ED 


here our number of Trotter slices is related to the number 
of qubits through T = 4 n. In reporting QMC times, we 
will report the number of sweeps each simulation took 
while going through the critical s region. There are n ■ T 
Metropolis steps per sweep, so the actual run-time of 
the algorithm depends polynomially on the number of 
sweeps. 

In Fig. [3 we show the full results of our QMC simu¬ 
lations, comparing the run-times of these algorithms to 
the corresponding There is a strong correlation be¬ 
tween these two quantities, which at least indicates some 
relation. The following sections will breakdown this data 
by a value and analyze it independently. 


A. Barriers Proportional to n° 5 

To start, we will focus on a = 0.5. Based on Fig. [2 
this size of barrier has QAO run-times that scale super- 
polynomially with n. Practically, we are able to run QMC 
simulations with n ranging up to ~ 220 qubits. For this 
regime of n, small n effects mask the superpolynomial 
scaling of the gap for c = 3 but not for c = 2. 

Note that c = 2 leads to smaller spectral gaps than 
c = 3 at fixed n. From trial and error, we found that the 
smaller gap sizes mean that the Trotter approximation 
needs to be better in order to get sensible results. Thus, 
for c = 2, T = 16n rather than the usual T = 4n. This 
necessity to improve the QMC for simulations with smaller 
gap sizes lends significant credence to the idea that the 
QMC algorithm depends heavily on the spectral gap itself. 

The QMC run-time (averaged over multiple simula¬ 
tions) as a function of g~ in is shown in Fig. [8] Notice 
that the data in this figure does not lie along a straight 
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l/(Minimum Gap Size) 2 

FIG. 8: QMC Sweeps vs. g m in for barrier size cn 0 ' 5 : The 
number of sweeps is increasing faster than a power law with 
the inverse gap size, indicating the our specific QMC algorithm 
is worse than QAO in this case. For c = 2 (green), n ranges 
from 84 to 172, and for c = 3 (blue), n ranges from 88 to 216. 



FIG. 9: QMC Sweeps vs. g m i n for barrier size cn 0 ' 4 : There 
appears to be a linear relationship here, indicating that QMC 
performance and QAO performance are polynomially related 
in this region. For c = 1 (red), n ranges from 184 to 320, for 
c = 2 (green), n ranges from 132 to 320, and for c = 3 (blue), 
n ranges from 116 to 224. 


line, so the QMC run-times seem to be increasing at a rate 
faster than polynomially in the inverse gap. This lack of 
a power law could be caused by three possible effects. 

It is possible that this means the QMC algorithm does 
indeed scale superpolynomially with c/~ in . An alterna¬ 
tive is that this curvature is due to small n effects that 
are still prevalent even for n in the several hundreds. Es¬ 
pecially for c = 3 and lower n, there is overlap between 
the initial s = 0 ground state distribution and the bar¬ 
rier, which could account for the apparent deviation from 
a power law here. Additionally, this curvature could be 
an indication of deficiencies in our QMC implementation 
specifically. As will be discussed in the next section, our 
algorithm has some notable approximations and simpli¬ 
fications that could be leading to this discrepancy. 


B. Barriers Proportional to n° 4 

For a = 0.4, the QMC simulations are able to go up 
to ~ 320 qubits. In this regime of n, small n effects 
mean that the gap is not superpolynomial for c = 3, 2 
(see Fig. GO) but it is for c = 1 (see Fig. 0. In Fig. M we 
have compared the QMC run-times directly to the spectral 
gap. Notice that in this case, there does seem to be a 
linear relationship between the log-log data. Many of 
the deficiencies in our specific implementation are less 
pronounced in this case than in the a = 0.5 case since the 
barrier is smaller. There is less overlap between the initial 
ground state and the barrier, which could also mean these 
simulations suffer less from small n effects than the a = 
0.5 simulations. 
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FIG. 10: QMC Sweeps vs. g m in for barrier size cn 0 ' 3 : There 
appears to be a linear relationship here, indicating that QMC 
performance and QAO performance are polynomially related 
in this region. For c = 1 (red), n ranges from 104 to 396, for 
c = 2 (green), n ranges from 104 to 660, and for c = 3 (blue), 
n ranges from 104 to 396. 

C. Barriers Proportional to n° 3 

Finally for a = 0.3, numerical diagonalization indi¬ 
cates the gap decreases polynomially in n for low and 
high n, no matter what c is chosen. Since the width of 
the barrier does not increase often for such a low scaling 
power a, the number of n accessible to the QMC sim¬ 
ulations is low here. Our data is displayed in Fig. [TO] 
Notice that there does seem to be a linear relationship 
on the log-log scale between inverse gap size and run-time 
here, though it is partially masked by the dearth of data 














FIG. 11: QMC Sweeps vs. s for barrier size 3n 0 ’ 3 at n = 116: 
This is averaged over 30 simulations. Notice that unlike Fig. [6] 
there is no noticeable spike here corresponding to tunneling. 

points. However, this does seem to indicate a polynomial 
relationship between QMC run-time and </“ in . 

Additionally, a plot of run-time versus s for higher 
powers, such as in Fig. [(3 shows a noticeable spike right at 
the tunneling location. For lower powers, such as a = 0.3 
as shown in Fig. 1111 there is no noticeable tunneling spike 
in the run-time. From our simulation results, it seems 
that the distinction between spikes and no spikes corre¬ 
sponds with the superpolynomial scaling cutoff we saw 
in the spectral gap in Section ITTT1 

VI. CONCLUSION 

First, in Section cm we numerically verified a folklore 
result [14jj about the relationship between n and the min¬ 
imum gap g m in- We showed that g m ; n scales polynomi- 
ally with n for barriers whose height and width grow like 
a < | but that for a > the minimum gap decreases 
faster than a power law. This indicates that QAO can 
succeed in finding the true ground state in polynomial 
time only for a < |. 

Our numerical results with Quantum Monte Carlo sim¬ 
ulations show that above a = there is a clear slowdown 
in the QMC algorithm (see Fig. [6]) whose location in s 
corresponds well with the location of the minimum gap 


in QAO. This slowdown all but disappears for lower a 
(see Fig. fllTl where the QMC algorithm has little trouble 
tunneling through the potential barrier. This is strong 
evidence that there is a correlation between spectral gap 
and QMC performance. 

Furthermore, in Section 0 we showed that there is in¬ 
deed a correlation between gap size and QMC run-time. 
For as less than we see data consistent with a polyno¬ 
mial relationship between QMC run-time and g~ in - This 
relationship is more difficult to discern for a > ^ with 
there seeming to be either a polynomial or superpolyno¬ 
mial relationship. The lack of a solid polynomial rela¬ 
tionship could be due to small n effects which are more 
prevalent in our simulations for higher a, or it could also 
be due to inadequacies in our QMC implementation rather 
than QMC algorithms in general 

Most notably our algorithm keeps a fixed As through¬ 
out its annealing schedule and relies on spending more 
time on each s value rather than decreasing the size of 
the s step. A more advanced algorithm could also dy¬ 
namically update s to move more slowly through problem 
regions. 

For the most part, our simulations also keep the num¬ 
ber of Trotter time steps T = An. While T = An is suffi¬ 
cient for the region of parameter space discussed in this 
article, it is possible that other Trotterization divisions 
would be more efficient 

Of course our work can also be extended by consid¬ 
ering different regions in parameter space of the Hamil¬ 
tonian. The scaling of the height and width are varied 
together using a in our analysis, but they can be varied 
independently. Additionally, the shape of the barrier can 
be made more complicated than the simple step used 
here. More generally, this procedure of applying QMC 
algorithms with annealing schedules can be used with 
other Hamiltonians to gain insight into the relationship 
between QAO and classical computing. 
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Appendix A: Derivation of Partition Function 

We will start with Eq. 0 and our goal will be to derive 
Eq. [9] as well as an estimator for our quantum mechanical 
ground state energy. Our first step will involve inserting 
our exponential approximation scheme so that (we will 
start using hats on operators to avoid confusion) 


Z = 


lim 

T—> oo 

x(°) 


E 

,...,rA T -!) 


W[x {T) \e~T' Hd e~T' Ho \x {T+1) ) , 

_T—0 


(Al) 


where the sums go over each x^ £ {0,1}". 

Here H a and H d are the off-diagonal and diagonal parts 
of the Hamiltonian, given by 


H d = Y M) \x)(x\ 


xG{0,l}" 


Ho=Y 

(x,y) 


(1 - s) 


\x)(y\- 


The sum in H a is over nearest neighbor sites (i.e. bit 
strings x and y that differ by one bit flip). Since H d is 
diagonal in the computational basis, we can just have 
it act on our basis states pulling out the eigenvalues 
H d {x ) = (l-s)f Ts/(|x|). 


Z = lim 

T — y oo 
T-l 


"T—1 

e n< 

A 0 ),...,®^- 1 ) Lr=o 


o -§,H d (x^) 


(A2) 


]V T) |e-^|z( r + 1 )} 


,T=0 


Next, we will claim that there is an orthonormal basis 
|that is the eigenbasis for H 0 , whose eigenvalues are 
H 0 (k < ' T )). We can insert a complete set of these states at 


ever time slice to get 


E E 


(A3) 


Z = lim 
T —too 

J^ 1 e -#^(,M) e -# ffo (W ) ) <;c ( T ) |fc (.) )(fcM | a .(r+l) ) 

. T—0 

To find these |fc) states, we just need to diagonalize 
H 0 . This operator can be represented by a translationally 
invariant matrix on an n dimensional hypercubic lattice 
(where each dimension is two sites long) with periodic 
boundary conditions and nearest neighbor interactions. 
These properties mean that the eigenstates of H 0 are 
simply the Brillouin Zone lattice sites. If we represent 
each Brillouin Zone lattice site using k £ {0,1}", then 
these lattice sites can be represented in the \x) basis by 


I k)= Y e 

xG{0,l} 


iirk-x 


x). 


(A4) 


Using standard Brillouin Zone methods for translation- 
ally invariant matrices, we can work out that the eigen¬ 
values of our off-diagonal Hamiltonian are 


H 0 (k) = - { ±-^Y{l-2k d ). 


(A5) 


d=1 


Furthermore, the overlap between \x) and |fc) states is 
given by 

(x\k) = (-1 (A6) 

Inserting Eqs. lAfil and lA5l back into our partition func¬ 
tion gives us 

TT-l 

o -e Hd (x^) 


Z = lim 

T— KX) 
T-l 


E 




(T-l) 


n e ~~ i 

,7=0 


(A7) 


n e n 

,T —0 fc(-r) d— 1 

We can rewrite TlLi YTd=i E fe w =0>r Focusing 
on just the important part and dropping the r labels in 
favor of labeling the two bit strings by x and y, we get 


H t >i 




n e 

d= 1 k d = 0,1 


T (1 2 S) (1 —2fcd) ^_^k d (x d -y d ) 


(AS) 


n[« ffii 


s) R (1-3) 


d—l 


Note that we have now eliminated the k variables entirely. 
Inserting this simplification lets us exactly recover Eq. [9] 


Z = lim 

T—»■ oo 


E 


A°). 


C (T-1) Lt=0 


T-l 


Y[ e -#((i-«)9+»/(l* w D) (A9) 


n 

n(^+(-!)■ 


(t) „(r+l) 


(1-3) 
g T 2 


d—l 
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Appendix B: Derivation of Energy Estimators 


Using Eq. IB31 the average becomes 


Next, we need to look at what the expectation value of 
a quantum operator is in the Trotter expanded formal¬ 
ism. By definition, we have 

(6) = i Tr{de-^}. (Bl) 



X 


lim 

T —>oo 


E K{* (r) }) 


(x^ T ~ 1 ^\e~T'^ d e~T'^ o 0\x^) 
(*( r -D \e~T^ d e~T’^° |x(°)} 


(B4) 


When we do the Trotter expansion we do not and 
should not expand O as we do the exponential. In fact af¬ 
ter the Trotter expansion, we will still only have one copy 
of O still, so the original copy of O will just be with one 
of the time slices. For convenience, we will put it with 
the very first time slice, so that after Trotterization, we 
are looking at: 


The QMC method will specifically use the average en¬ 
ergy: ^ H^ + (h o ). Starting with the 

operator is already acting on its eigenstates, so the aver¬ 
age becomes 




= lim — 
T—too Z 


E <■ 

t (»).lIT- 1 ) 




6|x (0) ) 


T—2 


-\{x( T) \e~T Hd e~T H °\x {T+l) ) 


r —0 


= lim — 

T—toc Z 


E 


- 1 ) 


(x^ T B|e tH<>q |2;( 0 )) 

(x( T - 1 }\e-T’H d e~'TH o \x(o)) 


T-l 


^[(z^le T H d e T Ho | x ( r + 1 )) 


T—0 


(B2) 



lim 

T—^oo 


c(°). 


E [ffd(* (0) )p({z (T) })] ■ (B5) 




In actual simulations, the estimator Hd{x^) —> 

*Er=o H d (x^) is used so that information from the 
entire time dimension can enter the statistics. 

Moving onto and focusing on just the relevant 

piece we have (replacing x^ T ^ 1 ' 1 —> x and x ^ — > y for 
notational convenience): 


Next consider the probability of obtaining a specific 
configuration, {x^}, of our n x T lattice of bits: 


p ({ lW }) s i 


T-l 

IK 

Lt— o 


"( r ) \p~ ~T^d ^H° I (i 


■+ 1 )\ 


(B3) 


{x\e~T , ^ d e~x fl ° H a \y) 
(x\e~T’^ d e~T’^° | y) 


we can insert k resolutions of the the identity in the top 
and bottom to get 


J 


Ea £ { o,i}" e ^ Ho{k) H 0 {k)(x\k)(k\y} 
Efc'e(o,i} n e~^ Ho ^ k,) (x\k')(k'\y) 


(i-s) 

2 


E 


6 t ( 2 5 Sdfcil 1 2fe <d Ep = i(l — 2k p )(— l) k ( x v "> 


feE{0,l} r 


Efc'e { o,i}" ^-^(-l) k '^-y) 

I- 


(B6) 


Next, we pull out what we can and switch 

EfcE{o,i} n EUi ~t Ild=i Efc d =o,i : 

In a given p element, the term in the product will be 
the same in the numerator and denominator if d ^ p, so 
the terms in the product cancel except in the case where 
d = p: 


(1 - s) -A eT^ - (_i)Up-yp) e - 

9 2^ 


f) (l-s) 

T 2 


2 ^ + (—1 )( x p Vp)e 


p d-g) 


(B7) 


Inserting Ea. lB7l into the off-diagonal energy estimator 


gives 


H n ) = lim 
T—>oo 


E ># w }) 


(B8) 




(1 — s) eT 2 — (— 1)^P X P >P 


E- 

p =i 


(0)_^.(T-1)^ 0 (1 ~s) 

e~T -2- 


er ' 


(1 -s) 




Again, we typically average over the result for the dif¬ 
ferent time slices in the actual simulation. 




























